knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE)

library(tidyverse)
library(skimr)
library(lubridate)
library(stringr)
library(nortest)
library(moments)
library(rstatix)
library(effectsize)
library(igraph)
library(tidygraph)
library(ggraph)
library(ggiraph)
library(broom)
library(robustbase)   
library(lmtest)       
library(performance)
library(see)
library(modelsummary)


# Formázott p-érték függvény
format_p <- function(p) {
  ifelse(p < 0.001, "< 0.001", as.character(round(p, 3)))}

Bevezetés

A gyógyszerek forgalomba hozatali engedélyezése összetett, több lépcsős szabályozási folyamat, amelyet az engedélyezést követően is folyamatos felügyelet és a dokumentáció ismtételt aktualizálása kísér. A forgalomba hozatali engedélyhez tartozó dokumentáció módosításai, az úgynevezett revíziók, a gyógyszer teljes életciklusa során tükrözhetik többek között a biztonsági jelzések, a klinikai bizonyítékok bővülésének, illetve a szabályozói eljárások sajátosságainak hatását. Emellett az engedélyezések mennyisége és összetétele (például terápiás területek szerint) időben is változhat, ami az gyógyszerfejlesztési trendek és a gyógyszercégek fókuszának változását mutathatják.

Az Európai Gyógyszerügynökség (European Medicines Agency, EMA) központi szerepet játszik a gyógyszerek engedélyezésében az Európai Unióban. Jelen elemzés célja az EMA gyógyszerengedélyezési adatbázisának részletes vizsgálata, amely négy fő kutatási kérdés köré szerveződik: 1. Az engedélyezések időbeli trendjeinek feltárása 1995 és 2023 között. 2. Az ATC-csoportok (Anatomical Therapeutic Chemical) közötti különbségek vizsgálata kulcsváltozók mentén 3. A revíziók számát befolyásoló tényezők azonosítása többváltozós modellezéssel 4. A központi idegrendszeri (ATC: N csoport) készítményeknél a terápiás területek és hatóanyagok kapcsolati szerkezetének hálózati szemléletű bemutatása.

Az elemzéshez az EMA nyilvánosan elérhető adatbázisát használtam, amely a TidyTuesday projekt keretében került publikálásra (2023. március 14.).

Adatok előkészítése

Az elemzés reprodukálható módon, webes forrásból betöltött adatbázison alapul, közvetlenül a TidyTuesday GitHub tárhelyéről került beolvasásra.

# Adat beolvasása a GitHub-ról
drugs <- 
  readr::read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/main/data/2023/2023-03-14/drugs.csv')
# Adatok áttekintése
dim(drugs)
## [1] 1988   28
names(drugs)
##  [1] "category"                                   
##  [2] "medicine_name"                              
##  [3] "therapeutic_area"                           
##  [4] "common_name"                                
##  [5] "active_substance"                           
##  [6] "product_number"                             
##  [7] "patient_safety"                             
##  [8] "authorisation_status"                       
##  [9] "atc_code"                                   
## [10] "additional_monitoring"                      
## [11] "generic"                                    
## [12] "biosimilar"                                 
## [13] "conditional_approval"                       
## [14] "exceptional_circumstances"                  
## [15] "accelerated_assessment"                     
## [16] "orphan_medicine"                            
## [17] "marketing_authorisation_date"               
## [18] "date_of_refusal_of_marketing_authorisation" 
## [19] "marketing_authorisation_holder_company_name"
## [20] "pharmacotherapeutic_group"                  
## [21] "date_of_opinion"                            
## [22] "decision_date"                              
## [23] "revision_number"                            
## [24] "condition_indication"                       
## [25] "species"                                    
## [26] "first_published"                            
## [27] "revision_date"                              
## [28] "url"
glimpse(drugs)
## Rows: 1,988
## Columns: 28
## $ category                                    <chr> "human", "human", "human",…
## $ medicine_name                               <chr> "Adcetris", "Nityr", "Ebva…
## $ therapeutic_area                            <chr> "Lymphoma, Non-Hodgkin;  H…
## $ common_name                                 <chr> "brentuximab vedotin", "ni…
## $ active_substance                            <chr> "brentuximab vedotin", "ni…
## $ product_number                              <chr> "002455", "004582", "00457…
## $ patient_safety                              <lgl> FALSE, FALSE, FALSE, FALSE…
## $ authorisation_status                        <chr> "authorised", "authorised"…
## $ atc_code                                    <chr> "L01XC12", "A16AX04", NA, …
## $ additional_monitoring                       <lgl> FALSE, FALSE, TRUE, TRUE, …
## $ generic                                     <lgl> FALSE, TRUE, FALSE, FALSE,…
## $ biosimilar                                  <lgl> FALSE, FALSE, FALSE, FALSE…
## $ conditional_approval                        <lgl> FALSE, FALSE, FALSE, FALSE…
## $ exceptional_circumstances                   <lgl> FALSE, FALSE, TRUE, FALSE,…
## $ accelerated_assessment                      <lgl> FALSE, FALSE, FALSE, FALSE…
## $ orphan_medicine                             <lgl> TRUE, FALSE, TRUE, FALSE, …
## $ marketing_authorisation_date                <date> 2012-10-25, 2018-07-26, 2…
## $ date_of_refusal_of_marketing_authorisation  <date> NA, NA, NA, NA, NA, NA, N…
## $ marketing_authorisation_holder_company_name <chr> "Takeda Pharma A/S", "Cycl…
## $ pharmacotherapeutic_group                   <chr> "Antineoplastic agents", "…
## $ date_of_opinion                             <date> 2012-07-19, 2018-05-31, 2…
## $ decision_date                               <date> 2022-11-17, 2023-03-10, 2…
## $ revision_number                             <dbl> 34, 4, 2, 3, 30, 24, 4, 18…
## $ condition_indication                        <chr> "Hodgkin lymphomaAdcetris …
## $ species                                     <chr> NA, NA, NA, NA, NA, NA, NA…
## $ first_published                             <dttm> 2018-07-25 13:58:00, 2018…
## $ revision_date                               <dttm> 2023-03-13 11:52:00, 2023…
## $ url                                         <chr> "https://www.ema.europa.eu…
skimr::skim(drugs)
Data summary
Name drugs
Number of rows 1988
Number of columns 28
_______________________
Column type frequency:
character 13
Date 4
logical 8
numeric 1
POSIXct 2
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
category 0 1.00 5 10 0 2 0
medicine_name 0 1.00 3 125 0 1976 0
therapeutic_area 285 0.86 4 400 0 669 0
common_name 4 1.00 4 220 0 1261 0
active_substance 1 1.00 4 823 0 1345 0
product_number 0 1.00 6 6 0 1932 0
authorisation_status 1 1.00 7 10 0 3 0
atc_code 28 0.99 3 18 0 1074 0
marketing_authorisation_holder_company_name 4 1.00 4 65 0 615 0
pharmacotherapeutic_group 34 0.98 7 174 0 365 0
condition_indication 12 0.99 18 7597 0 1886 0
species 1709 0.14 4 67 0 59 0
url 0 1.00 53 148 0 1988 0

Variable type: Date

skim_variable n_missing complete_rate min max median n_unique
marketing_authorisation_date 60 0.97 1995-10-20 2023-02-20 2013-06-09 1127
date_of_refusal_of_marketing_authorisation 1913 0.04 2004-09-07 2022-04-29 2013-04-25 67
date_of_opinion 779 0.61 1995-07-12 2022-12-15 2016-07-21 389
decision_date 45 0.98 1998-08-20 2023-03-10 2022-02-16 815

Variable type: logical

skim_variable n_missing complete_rate mean count
patient_safety 0 1 0.01 FAL: 1977, TRU: 11
additional_monitoring 0 1 0.19 FAL: 1601, TRU: 387
generic 0 1 0.16 FAL: 1673, TRU: 315
biosimilar 0 1 0.05 FAL: 1896, TRU: 92
conditional_approval 0 1 0.02 FAL: 1940, TRU: 48
exceptional_circumstances 0 1 0.02 FAL: 1940, TRU: 48
accelerated_assessment 0 1 0.02 FAL: 1940, TRU: 48
orphan_medicine 0 1 0.08 FAL: 1826, TRU: 162

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
revision_number 96 0.95 13.53 11.65 0 4.75 11 19 89 ▇▃▁▁▁

Variable type: POSIXct

skim_variable n_missing complete_rate min max median n_unique
first_published 0 1.00 1998-08-20 00:00:00 2023-03-09 18:50:00 2018-04-14 23:14:30 1760
revision_date 29 0.99 2000-07-17 02:00:00 2023-03-13 11:52:00 2022-05-11 11:58:00 1932
# Folytonos változók kialakítása
## Gyógyszer piacon töltött idejének (days_on_market változó) kiszámítása napokban (Az adatbázis feltöltési dátuma: 2023.03.13.)
drugs <- drugs %>%
  mutate(
    days_on_market = as.numeric(ymd("2023-03-13") - ymd(marketing_authorisation_date)))

## Ellenörzés - days_on_market
summary(drugs$days_on_market)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      21    1818    3564    3948    5655   10006      60
## Eloszlás - days_on_market
ggplot(drugs, aes(x = days_on_market)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  labs(x = "Napok száma a piacon", y = "Gyakoriság") +
  theme_minimal()

A „piacon töltött időt” (days_on_market változó) a forgalomba hozatali engedély dátuma és az adatbázis feltöltési dátuma (2023-03-13) közötti napok száma alapján képeztem, ezzel biztosítva, hogy a mutató ne függjön a futtatás aktuális dátumától.

Adatok feldolgozása

Minta

Az adatbázis összesen 1988 gyógyszert tartalmaz, amelyek közül 1706 (85.8%) humán és 282 (14.2%) állatgyógyászati készítmény. Az adatbázis 28 változót tartalmaz, beleértve a gyógyszer nevét, terápiás területét, hatóanyagát, ATC kódját, engedélyezési dátumát és státuszát, valamint olyan különböző jellemzőket a gyógyszerekkel kapcsolatban, mint pl. generikus, orphan, gyorsított eljárás stb.

A gyógyszerek engedélyezési dátuma 1995. október 20. és 2023. február 20. között oszlik meg (medián: 2013. június 9.). Az engedélyezési státusz tekintetében a gyógyszerek 79.1%-a (n = 1573) rendelkezik érvényes engedéllyel, 18.0% (n = 357) visszavont, és 2.9% (n = 57) elutasított státuszú.

A speciális státuszok tekintetében a gyógyszerek 15.8%-a (n = 315) generikus, 8.2%-a (n = 162) orphan gyógyszer és 2.4%-a (n = 48) részesült gyorsított eljárásban.

A gyógyszerek piacon töltött idejének eloszlása a teljes mintában széles tartományt fed le (Min: 21 nap; Mdn: 3564 nap; M: 3948 nap; Max: 10006 nap; NA: 60).

# Kategórikus változók gyakorisága és százalékos eloszlása
vars <- c(
  "category", "patient_safety", "authorisation_status", "generic",
  "conditional_approval", "exceptional_circumstances",
  "accelerated_assessment", "orphan_medicine")

cat_tables <- drugs %>%
  select(all_of(vars)) %>%
  lapply(function(x) {
    tab <- table(x, useNA = "ifany")
    pct <- prop.table(tab) * 100
    data.frame(
      value   = names(tab),
      n       = as.integer(tab),
      percent = round(as.numeric(pct), 1),
      row.names = NULL)})

cat_tables$category
##        value    n percent
## 1      human 1706    85.8
## 2 veterinary  282    14.2
cat_tables$patient_safety
##   value    n percent
## 1 FALSE 1977    99.4
## 2  TRUE   11     0.6
cat_tables$authorisation_status
##        value    n percent
## 1 authorised 1573    79.1
## 2    refused   57     2.9
## 3  withdrawn  357    18.0
## 4       <NA>    1     0.1
cat_tables$generic
##   value    n percent
## 1 FALSE 1673    84.2
## 2  TRUE  315    15.8
cat_tables$conditional_approval
##   value    n percent
## 1 FALSE 1940    97.6
## 2  TRUE   48     2.4
cat_tables$exceptional_circumstances
##   value    n percent
## 1 FALSE 1940    97.6
## 2  TRUE   48     2.4
cat_tables$accelerated_assessment
##   value    n percent
## 1 FALSE 1940    97.6
## 2  TRUE   48     2.4
cat_tables$orphan_medicine
##   value    n percent
## 1 FALSE 1826    91.9
## 2  TRUE  162     8.1
# Piacon töltött idő (days_on_market) leíró statisztikája
summary(drugs$days_on_market)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      21    1818    3564    3948    5655   10006      60
# Gyógyszerek engedélyezésének száma évente
drugs %>%
  filter(!is.na(marketing_authorisation_date)) %>%
  mutate(year = year(marketing_authorisation_date)) %>%
  count(year) %>%
  ggplot(aes(year, n)) +
  geom_line() +
  scale_x_continuous(breaks = seq(1995, 2022, by = 1), minor_breaks = NULL, limits = c(NA, 2022)) +
  labs(title = "Gyógyszer-engedélyezések száma évente", x = "Év", y = "n") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Top 10 gyártó
drugs %>%
  count(marketing_authorisation_holder_company_name, sort = TRUE) %>%
  head(10)
## # A tibble: 10 × 2
##    marketing_authorisation_holder_company_name     n
##    <chr>                                       <int>
##  1 Accord Healthcare S.L.U.                       58
##  2 Novartis Europharm Limited                     58
##  3 Pfizer Europe MA EEIG                          43
##  4 Zoetis Belgium SA                              40
##  5 AstraZeneca AB                                 36
##  6 Boehringer Ingelheim Vetmedica GmbH            36
##  7 Merck Sharp & Dohme B.V.                       32
##  8 Teva B.V.                                      32
##  9 Intervet International BV                      31
## 10 Eli Lilly Nederland B.V.                       30

Gyógyszer-engedélyezési trendek

Az 1. ábra a gyógyszer-engedélyezések számának alakulását mutatja az 1995-2022 közötti időszakban. A trendvonal jelentős fluktuációt mutat, kiugró csúcsokkal 2001 (n ≈ 90), 2017 (n ≈ 110) és 2021 körül (n ≈ 120). A COVID-19 pandémia időszakában (2020-2021) a gyógyszer-engedélyezések száma jelentősen megugrott, ami feltételezhetően a vakcinák és vírusellenes szerek gyorsított engedélyezésének tudható be elsősorban.

A forgalmazási engedély jogosultjai közül a legtöbb gyógyszert az Accord Healthcare S.L.U. és a Novartis Europharm Limited jegyzi (egyaránt n = 58), amelyeket a Pfizer Europe MA EEIG (n = 43) és a Zoetis Belgium SA (n = 40) követ.

ATC-rendszer

A humán készítmények ATC-főcsoport szerinti bontása lehetővé teszi annak áttekintését, hogy az engedélyezések elsősorban mely terápiás területekhez kapcsolódtak, illetve hogyan alakultak ezen időszakban. Az ATC-főcsoportok szerinti évekre vontott ábra a főcsoportok közötti nagyságrendi különbségeket és az időbeli dinamikát együttesen teszi láthatóvá, így a későbbi, célzottabb összehasonlítások megalapozására szolgálhat.

# Human gyógyszerek ATC kód csoportonként
drugs %>%
  filter(category == "human", !is.na(marketing_authorisation_date), !is.na(atc_code)) %>%
  mutate(
    year = year(marketing_authorisation_date),
    atc_group = str_extract(atc_code, "^[A-Z]")) %>%
  count(year, atc_group) %>%
  ggplot(aes(year, n)) +
  geom_line() +
  facet_wrap(~atc_group, scales = "free_x") +
  scale_x_continuous(breaks = seq(1995, 2022, by = 5), limits = c(1995, 2022), minor_breaks = NULL) +
  labs(
    title = "Humán gyógyszer-engedélyezések száma évente ATC-rendszer szerint",
    x = "Év",
    y = "n",
    caption = "Note. ATC besorolás: A - Tápcsatorna és anyagcsere; B - Vér és vérképzőszervek; C - Cardiovascularis rendszer; D - Bőrgyógyászati készítmények; G - Urogenitalis rendszer és nemi hormonok;\nH - Szisztémás hormonkészítmények; J - Szisztémás fertőzésellenes szerek; L - Daganatellenes és immunmoduláns szerek; M - Váz- és izomrendszer; N - Központi idegrendszer;\nP - Parazitaellenes készítmények; R - Légzőrendszer; S - Érzékszervek; V - Egyéb") +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.text = element_text(size = 10),
    plot.caption = element_text(hjust = 0, size = 9))

# Heatmap - ATC csoportok évenete
drugs %>%
  filter(category == "human", !is.na(marketing_authorisation_date), !is.na(atc_code)) %>%
  mutate(
    year = year(marketing_authorisation_date),
    atc_group = str_extract(atc_code, "^[A-Z]")) %>%
  filter(year <= 2022) %>%
  count(year, atc_group) %>%
  group_by(year) %>%
  mutate(percent = n / sum(n) * 100) %>%
  ungroup() %>%
  ggplot(aes(x = atc_group, y = factor(year), fill = percent)) +
  geom_tile(color = "white") +
  geom_text(aes(label = ifelse(n > 0, n, "")), size = 2.5) +
  scale_fill_gradient(low = "white", high = "steelblue", na.value = "white") +
  labs(
    title = "Humán gyógyszer-engedélyezésének száma ATC-rendszer szerint",
    x = "ATC csoport",
    y = "Év",
    fill = "Arány (%)",
    caption = "Note. ATC besorolás: A - Tápcsatorna és anyagcsere; B - Vér és vérképzőszervek; C - Cardiovascularis rendszer; D - Bőrgyógyászati készítmények; G - Urogenitalis rendszer és nemi hormonok; \nH - Szisztémás hormonkészítmények; J - Szisztémás fertőzésellenes szerek; L - Daganatellenes és immunmoduláns szerek; M - Váz- és izomrendszer; N - Központi idegrendszer; \nP - Parazitaellenes készítmények; R - Légzőrendszer; S - Érzékszervek; V - Egyéb") +
  theme_minimal() +
  theme(
    plot.caption = element_text(hjust = 0, size = 9),
    panel.grid = element_blank())

# ATC csoportok gyakorisága a teljes időszakban
drugs %>%
  filter(category == "human", !is.na(atc_code)) %>%
  mutate(atc_group = str_extract(atc_code, "^[A-Z]")) %>%
  count(atc_group, sort = TRUE) %>%
  mutate(percent = round(n / sum(n) * 100, 1))
## # A tibble: 14 × 3
##    atc_group     n percent
##    <chr>     <int>   <dbl>
##  1 L           463    27.5
##  2 J           253    15.1
##  3 A           201    12  
##  4 N           186    11.1
##  5 B           135     8  
##  6 C            95     5.7
##  7 V            72     4.3
##  8 R            69     4.1
##  9 G            59     3.5
## 10 M            54     3.2
## 11 H            40     2.4
## 12 S            34     2  
## 13 D            18     1.1
## 14 P             2     0.1

Az ATC-rendszer szerinti megoszlás

Az ATC-rendszer szerinti elemzés a humán gyógyszerekre korlátozódott. A 14 fő ATC csoportból a leggyakoribb az L csoport (daganatellenes és immunmoduláns szerek; n = 463, 27.5%), amelyet a J csoport (szisztémás fertőzésellenes szerek; n = 253, 15.1%), az A csoport (tápcsatorna és anyagcsere; n = 201, 12.0%) és az N csoport (központi idegrendszer; n = 186, 11.1%) követ. Ez a négy csoport együttesen a humán gyógyszerek 65.7%-át teszi ki.

A heatmap-vizualizáció (2. ábra) az ATC csoportok éves megoszlását szemlélteti. Az L csoport dominanciája különösen az elmúlt évtizedben vált hangsúlyossá, ami feltehetően a precíziós onkológia és az immunterápiák fejlődésével magyarázható.

A további elemzéseket a négy leggyakoribb ATC csoportra (L, J, A, N) szűkítettem, amelyek együttesen 1103 gyógyszert foglalnak magukban.

Statisztikai elemzések - TOP4 ATC

# Top 4 ATC csoport kiszűrése
top4_atc <- drugs %>%
  filter(category == "human", !is.na(atc_code)) %>%
  mutate(atc_group = str_extract(atc_code, "^[A-Z]")) %>%
  filter(atc_group %in% c("L", "J", "A", "N"))

# Leíró statisztika - revision_number
top4_atc %>%
  group_by(atc_group) %>%
  summarise(
    n = n(),
    hianyzo = sum(is.na(revision_number)),
    atlag = mean(revision_number, na.rm = TRUE),
    sd = sd(revision_number, na.rm = TRUE),
    min = min(revision_number, na.rm = TRUE),
    q25 = quantile(revision_number, 0.25, na.rm = TRUE),
    median = median(revision_number, na.rm = TRUE),
    q75 = quantile(revision_number, 0.75, na.rm = TRUE),
    max = max(revision_number, na.rm = TRUE))
## # A tibble: 4 × 10
##   atc_group     n hianyzo atlag    sd   min   q25 median   q75   max
##   <chr>     <int>   <int> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1 A           201       7  14.1  9.51     0     6     14  21      37
## 2 J           253       6  17.7 14.2      0     7     14  26.5    68
## 3 L           463      31  14.5 12.9      0     5     11  20      89
## 4 N           186       8  15.0 11.6      0     7     13  21      58
# Leíró statisztika - days_on_market
top4_atc %>%
  group_by(atc_group) %>%
  summarise(
    n = n(),
    hianyzo = sum(is.na(days_on_market)),
    atlag = mean(days_on_market, na.rm = TRUE),
    sd = sd(days_on_market, na.rm = TRUE),
    min = min(days_on_market, na.rm = TRUE),
    q25 = quantile(days_on_market, 0.25, na.rm = TRUE),
    median = median(days_on_market, na.rm = TRUE),
    q75 = quantile(days_on_market, 0.75, na.rm = TRUE),
    max = max(days_on_market, na.rm = TRUE))
## # A tibble: 4 × 10
##   atc_group     n hianyzo atlag    sd   min   q25 median   q75   max
##   <chr>     <int>   <int> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1 A           201       4 4182. 2524.    94 2163   3979  5908   9813
## 2 J           253       2 4300. 2774.    98 2008.  4042  6208.  9805
## 3 L           463      20 3130. 2455.    26 1138   2482  4676   9968
## 4 N           186       8 3987. 2420.   126 2267   3808. 5258.  9772
# ATC csoport nevek hozzáadása
top4_atc <- top4_atc %>%
  mutate(atc_name = case_when(
    atc_group == "L" ~ "L - Daganatellenes és\nimmunomoduláns szerek",
    atc_group == "J" ~ "J - Szisztémás\nfertőzésellenes szerek",
    atc_group == "A" ~ "A - Tápcsatorna\nés anyagcsere",
    atc_group == "N" ~ "N - Központi\nidegrendszer"
  ))

# Violin plot - revision_number
ggplot(top4_atc, aes(x = atc_name, y = revision_number, fill = atc_name)) +
  geom_violin(alpha = 0.5) +
  geom_boxplot(width = 0.2, outlier.shape = NA) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "red") +
  labs(
    title = "Revíziók száma ATC csoportonként",
    x = "ATC csoport",
    y = "Revíziók száma"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

# Violin plot - days_on_market
ggplot(top4_atc, aes(x = atc_name, y = days_on_market, fill = atc_name)) +
  geom_violin(alpha = 0.5) +
  geom_boxplot(width = 0.2, outlier.shape = NA) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "red") +
  labs(
    title = "Piacon eltöltött napok száma ATC csoportonként",
    x = "ATC csoport",
    y = "Napok száma"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

Leíró statisztika

A revíziók számának leíró statisztikáját az 1. táblázat foglalja össze. A revíziók átlagos száma a négy ATC csoportban 14.1 és 17.7 között mozgott. A legmagasabb átlagot a J csoport mutatta (M = 17.7, SD = 14.2), míg a legalacsonyabbat az A csoport (M = 14.1, SD = 9.5).

A piacon töltött idő tekintetében az L csoport gyógyszerei átlagosan rövidebb ideje vannak a piacon (M = 3130 nap, SD = 2455), mint az A (M = 4182 nap, SD = 2524) és J (M = 4300 nap, SD = 2774) csoport készítményei.

A violin plotok (3. és 4. ábra) vizuálisan is szemléltetik a csoportok közötti különbségeket.

# Függő változó: revision_number
## Shapiro-Wilk teszt - revision_number
top4_atc %>%
  group_by(atc_name) %>%
  summarise(
    W = round(shapiro.test(revision_number)$statistic, 3),
    p_value = format_p(shapiro.test(revision_number)$p.value))
## # A tibble: 4 × 3
##   atc_name                                           W p_value
##   <chr>                                          <dbl> <chr>  
## 1 "A - Tápcsatorna\nés anyagcsere"               0.959 < 0.001
## 2 "J - Szisztémás\nfertőzésellenes szerek"       0.922 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.87  < 0.001
## 4 "N - Központi\nidegrendszer"                   0.919 < 0.001
## Lilliefors teszt - revision_number
top4_atc %>%
  group_by(atc_name) %>%
  summarise(
    D = round(lillie.test(revision_number)$statistic, 3),
    p_value = format_p(lillie.test(revision_number)$p.value))
## # A tibble: 4 × 3
##   atc_name                                           D p_value
##   <chr>                                          <dbl> <chr>  
## 1 "A - Tápcsatorna\nés anyagcsere"               0.108 < 0.001
## 2 "J - Szisztémás\nfertőzésellenes szerek"       0.138 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.132 < 0.001
## 4 "N - Központi\nidegrendszer"                   0.123 < 0.001
# Függő változó: days_on_market
## Shapiro-Wilk teszt - days_on_market
top4_atc %>%
  group_by(atc_name) %>%
  summarise(
    W = round(shapiro.test(days_on_market)$statistic, 3),
    p_value = format_p(shapiro.test(days_on_market)$p.value))
## # A tibble: 4 × 3
##   atc_name                                           W p_value
##   <chr>                                          <dbl> <chr>  
## 1 "A - Tápcsatorna\nés anyagcsere"               0.965 < 0.001
## 2 "J - Szisztémás\nfertőzésellenes szerek"       0.945 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.911 < 0.001
## 4 "N - Központi\nidegrendszer"                   0.958 < 0.001
## Lilliefors teszt - days_on_market
top4_atc %>%
  group_by(atc_name) %>%
  summarise(
    D = round(lillie.test(days_on_market)$statistic, 3),
    p_value = format_p(lillie.test(days_on_market)$p.value))
## # A tibble: 4 × 3
##   atc_name                                           D p_value
##   <chr>                                          <dbl> <chr>  
## 1 "A - Tápcsatorna\nés anyagcsere"               0.069 0.025  
## 2 "J - Szisztémás\nfertőzésellenes szerek"       0.094 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.116 < 0.001
## 4 "N - Központi\nidegrendszer"                   0.062 0.091
# Ferdeség és csúcsosság
top4_atc %>%
  group_by(atc_name) %>%
  summarise(
    skewness_revision = round(skewness(revision_number, na.rm = TRUE), 3),
    kurtosis_revision = round(kurtosis(revision_number, na.rm = TRUE), 3),
    skewness_days = round(skewness(days_on_market, na.rm = TRUE), 3),
    kurtosis_days = round(kurtosis(days_on_market, na.rm = TRUE), 3))
## # A tibble: 4 × 5
##   atc_name       skewness_revision kurtosis_revision skewness_days kurtosis_days
##   <chr>                      <dbl>             <dbl>         <dbl>         <dbl>
## 1 "A - Tápcsato…             0.357              2.19         0.286          2.12
## 2 "J - Szisztém…             0.93               3.34         0.378          2.02
## 3 "L - Daganate…             1.55               6.48         0.912          3.05
## 4 "N - Központi…             1.07               4.00         0.528          2.73
# Vizuális módszerek
## Hisztogram - revision_number
ggplot(top4_atc, aes(x = revision_number, fill = atc_name)) +
  geom_histogram(bins = 30, color = "white") +
  facet_wrap(~atc_name, scales = "free") +
  labs(
    title = "Revíziók számának eloszlása ATC csoportonként",
    x = "Revíziók száma",
    y = "Gyakoriság") +
  theme_minimal() +
  theme(legend.position = "none")

## Q-Q plot - revision_number
ggplot(top4_atc, aes(sample = revision_number, color = atc_name)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~atc_name, scales = "free") +
  labs(
    title = "Q-Q plot: Revíziók száma ATC csoportonként",
    x = "Elméleti kvantilisek",
    y = "Minta kvantilisek") +
  theme_minimal() +
  theme(legend.position = "none")

## Hisztogram - days_on_market
ggplot(top4_atc, aes(x = days_on_market, fill = atc_name)) +
  geom_histogram(bins = 30, color = "white") +
  facet_wrap(~atc_name, scales = "free") +
  labs(
    title = "Piacon eltöltött napok eloszlása ATC csoportonként",
    x = "Napok száma",
    y = "Gyakoriság") +
  theme_minimal() +
  theme(legend.position = "none")

## Q-Q plot - days_on_market
ggplot(top4_atc, aes(sample = days_on_market, color = atc_name)) +
  stat_qq() +
  stat_qq_line() +
  facet_wrap(~atc_name, scales = "free") +
  labs(
    title = "Q-Q plot: Piacon eltöltött napok ATC csoportonként",
    x = "Elméleti kvantilisek",
    y = "Minta kvantilisek") +
  theme_minimal() +
  theme(legend.position = "none")

A normalitás vizsgálata

A paraméteres statisztikai eljárások alkalmazhatóságának előfeltételeként a függő változók normalitását Shapiro-Wilk és Lilliefors (Kolmogorov-Smirnov korrekció) tesztekkel vizsgáltam.

A revíziók számának eloszlása a Shapiro-Wilk teszt alapján egyik ATC csoportban sem követte a normális eloszlást (p < .001). A Lilliefors teszt eredményei megerősítették ezt a mintázatot (p < .001). A ferdeségi mutatók pozitív értékeket vettek fel (skewness = 0.36–1.55), ami jobbra ferde eloszlásokra utal. A csúcsossági értékek (kurtosis = 2.19–6.48) az L csoportban kifejezetten leptokurtikus eloszlást jeleztek.

A piacon töltött idő esetében hasonló eredmények lárhatók. A Shapiro-Wilk teszt valamennyi csoportban szignifikáns eltérést mutatott a normalitástól (p < .001). A Lilliefors teszt az N csoport kivételével (p = .091) szintén szignifikáns eredményeket adott.

A Q-Q plotok (5-8. ábra) vizuálisan is alátámasztják a normalitástól való eltérést, különösen az eloszlások szélein megfigyelhető szisztematikus eltérésekkel.

Tekintettel a normalitás sérülésére, a csoportok összehasonlításához nem-paraméteres statisztikai eljárásokat alkalmaztam.

# Függő változó: revison_number
## Kruskal-Wallis teszt - revision_number
top4_atc %>%
  kruskal_test(revision_number ~ atc_group) %>%
  mutate(p = format_p(p))
## # A tibble: 1 × 6
##   .y.                 n statistic    df p     method        
##   <chr>           <int>     <dbl> <int> <chr> <chr>         
## 1 revision_number  1103      9.39     3 0.025 Kruskal-Wallis
## Kruskal-Wallis hatásméret (eta²) - revision_number
top4_atc %>%
  kruskal_effsize(revision_number ~ atc_group, ci = TRUE)
## # A tibble: 1 × 7
##   .y.                 n effsize conf.low conf.high method  magnitude
## * <chr>           <int>   <dbl>    <dbl>     <dbl> <chr>   <ord>    
## 1 revision_number  1103 0.00581 -0.00073      0.02 eta2[H] small
## Dunn-teszt Holm korrekcióval - revision_number
dunn_revision <- dunn_test(revision_number ~ atc_group, data = top4_atc, p.adjust.method = "holm")
dunn_revision %>%
  mutate(p.adj = format_p(p.adj))
## # A tibble: 6 × 9
##   .y.             group1 group2    n1    n2 statistic       p p.adj p.adj.signif
##   <chr>           <chr>  <chr>  <int> <int>     <dbl>   <dbl> <chr> <chr>       
## 1 revision_number A      J        194   247    1.60   0.110   0.552 ns          
## 2 revision_number A      L        194   432   -1.05   0.292   0.833 ns          
## 3 revision_number A      N        194   178    0.0537 0.957   0.957 ns          
## 4 revision_number J      L        247   432   -3.06   0.00219 0.013 *           
## 5 revision_number J      N        247   178   -1.50   0.133   0.552 ns          
## 6 revision_number L      N        432   178    1.09   0.278   0.833 ns
## Cliff's delta függvény 95% CI-vel
calc_cliffs_delta <- function(data, var, group1, group2) {
  x <- data %>% filter(atc_group == group1) %>% pull({{var}})
  y <- data %>% filter(atc_group == group2) %>% pull({{var}})
  result <- effectsize::cliffs_delta(x, y, ci = 0.95)
  tibble(
    group1 = group1,
    group2 = group2,
    cliffs_delta = round(result$r, 3),
    ci_low = round(result$CI_low, 3),
    ci_high = round(result$CI_high, 3)
  )
}

## Cliff's delta páronként - revision_number
pairs <- combn(unique(top4_atc$atc_group), 2, simplify = FALSE)
map_df(pairs, ~calc_cliffs_delta(top4_atc, revision_number, .x[1], .x[2]))
## # A tibble: 6 × 5
##   group1 group2 cliffs_delta ci_low ci_high
##   <chr>  <chr>         <dbl>  <dbl>   <dbl>
## 1 L      A            -0.055 -0.152   0.043
## 2 L      J            -0.137 -0.224  -0.048
## 3 L      N            -0.059 -0.158   0.042
## 4 A      J            -0.094 -0.201   0.014
## 5 A      N            -0.001 -0.118   0.116
## 6 J      N             0.088 -0.023   0.197
# Függő változó: days_on_market
## Kruskal-Wallis teszt - days_on_market
top4_atc %>%
  kruskal_test(days_on_market ~ atc_group) %>%
  mutate(p = format_p(p))
## # A tibble: 1 × 6
##   .y.                n statistic    df p       method        
##   <chr>          <int>     <dbl> <int> <chr>   <chr>         
## 1 days_on_market  1103      48.2     3 < 0.001 Kruskal-Wallis
## Kruskal-Wallis hatásméret (eta²) - days_on_market
top4_atc %>%
  kruskal_effsize(days_on_market ~ atc_group, ci = TRUE)
## # A tibble: 1 × 7
##   .y.                n effsize conf.low conf.high method  magnitude
## * <chr>          <int>   <dbl>    <dbl>     <dbl> <chr>   <ord>    
## 1 days_on_market  1103  0.0412     0.02      0.07 eta2[H] small
## Dunn-teszt Holm korrekcióval - days_on_market
dunn_days <- dunn_test(days_on_market ~ atc_group, data = top4_atc, p.adjust.method = "holm")
dunn_days %>%
  mutate(p.adj = format_p(p.adj))
## # A tibble: 6 × 9
##   .y.            group1 group2    n1    n2 statistic        p p.adj p.adj.signif
##   <chr>          <chr>  <chr>  <int> <int>     <dbl>    <dbl> <chr> <chr>       
## 1 days_on_market A      J        197   251     0.105  9.17e-1 1     ns          
## 2 days_on_market A      L        197   443    -5.15   2.58e-7 < 0.… ****        
## 3 days_on_market A      N        197   178    -0.551  5.82e-1 1     ns          
## 4 days_on_market J      L        251   443    -5.71   1.13e-8 < 0.… ****        
## 5 days_on_market J      N        251   178    -0.683  4.94e-1 1     ns          
## 6 days_on_market L      N        443   178     4.33   1.50e-5 < 0.… ****
## Cliff's delta páronként - days_on_market
map_df(pairs, ~calc_cliffs_delta(top4_atc, days_on_market, .x[1], .x[2]))
## # A tibble: 6 × 5
##   group1 group2 cliffs_delta ci_low ci_high
##   <chr>  <chr>         <dbl>  <dbl>   <dbl>
## 1 L      A            -0.256 -0.344  -0.163
## 2 L      J            -0.254 -0.336  -0.169
## 3 L      N            -0.23  -0.323  -0.133
## 4 A      J            -0.013 -0.12    0.095
## 5 A      N             0.04  -0.077   0.156
## 6 J      N             0.047 -0.064   0.157

Csoportok összehasonlítása

Revíziók száma

A négy ATC csoport közötti különbséget Kruskal-Wallis teszttel vizsgáltam. Az eredmények szignifikáns különbséget mutattak a csoportok között (H(3) = 9.39, p = .025). A hatásméret ugyanakkor kicsi volt (η² = 0.006, 95% CI [-0.001, 0.020]). A 95% CI alsó határa 0 közelébe esik (a becslési eljárás miatt numerikusan enyhén negatív), ami a hatás elhanyagolható mértékével konzisztens. A Dunn-féle post-hoc teszt Holm-korrekcióval egyetlen szignifikáns páronkénti különbséget tárt fel. A J csoport (szisztémás fertőzésellenes szerek) szignifikánsan magasabb revíziószámmal rendelkezett, mint az L csoport (daganatellenes szerek), (padj = .013). A Cliff’s delta hatásméret mutató eredménye alapján kis hatásméret írható le a két csoport között (d = -0.137, 95% CI [-0.224, -0.048]).

Piacon töltött idő

A piacon töltött idő tekintetében a Kruskal-Wallis teszt erősen szignifikáns különbséget mutatott a csoportok között (H(3) = 48.2, p < .001). A hatásméret kis-közepes mértékű volt (η² = 0.041, 95% CI [0.020, 0.070]). A Dunn-féle post-hoc teszt három szignifikáns páronkénti különbséget tárt fel (padj < .001). Az L csoport gyógyszerei szignifikánsan rövidebb ideje vannak a piacon, mint az A csoport (d = -0.256), a J csoport (d = -0.254) és az N csoport (d = -0.230) készítményei. Ez a mintázat összhangban van az onkológiai területen tapasztalható intenzívebb kutatás-fejlesztési aktivitással és a gyorsabb gyógyszercserélődéssel.

Hálózatelemzés

Adatok előkészítése

# A therapeutic_area változó tördelése
drugs_clean <- drugs %>%
  separate_rows(therapeutic_area, sep = "\\s*;\\s*") %>% 
  mutate(therapeutic_area = trimws(therapeutic_area))

# N ATC csoport kiszűrése (központi idegrendszer)
drugs_N <- drugs_clean %>%
  filter(!is.na(atc_code)) %>%
  mutate(atc_group = str_extract(atc_code, "^[A-Z]")) %>%
  filter(atc_group == "N") %>%
  filter(!is.na(therapeutic_area))

# Összefoglaló statisztika
drugs_N %>%
  summarise(
    n_drugs = n_distinct(medicine_name),
    n_areas = n_distinct(therapeutic_area),
    n_substances = n_distinct(active_substance),
    n_companies = n_distinct(marketing_authorisation_holder_company_name))
## # A tibble: 1 × 4
##   n_drugs n_areas n_substances n_companies
##     <int>   <int>        <int>       <int>
## 1     185      70          101         111

A hálózatelemzés a központi idegrendszeri (N) ATC csoportra fókuszált, amely a pszichiátriai és neurológiai gyógyszereket foglalja magában. A vizsgálat célja a terápiás területek és hatóanyagok közötti kapcsolatok feltérképezése volt.

Az N csoportban összesen 185 egyedi gyógyszer, 70 terápiás terület, 101 hatóanyag és 111 forgalmazó cég volt azonosítható. A therapeutic_area változó pontosvesszővel elválasztott többszörös értékeket tartalmazott, amelyeket külön sorokra bontottam a hálózatépítéséhez.

Szűrés és hálózat építése

# Csak a gyakoribb terápiás területek és hatóanyagok (min. 3 előfordulás)
top_areas_N <- drugs_N %>%
  count(therapeutic_area) %>%
  filter(n >= 3) %>%
  pull(therapeutic_area)

top_substances_N <- drugs_N %>%
  count(active_substance) %>%
  filter(n >= 3) %>%
  pull(active_substance)

# Szűrt adatok
drugs_N_filtered <- drugs_N %>%
  filter(therapeutic_area %in% top_areas_N | active_substance %in% top_substances_N)

# Élsúlyok hozzáadása (hány gyógyszer kapcsolja össze őket)
edges_weighted <- drugs_N_filtered %>%
  count(therapeutic_area, active_substance, name = "weight") %>%
  rename(from = therapeutic_area, to = active_substance)

# ATC alcsoport kinyerése (N01, N02, stb.)
atc_subgroups <- drugs_N_filtered %>%
  select(therapeutic_area, active_substance, atc_code) %>%
  mutate(atc_subgroup = str_extract(atc_code, "^N\\d{2}")) %>%
  distinct(active_substance, atc_subgroup) %>%
  group_by(active_substance) %>%
  slice(1) %>%
  ungroup()

# Csúcsok létrehozása típussal és ATC alcsoporttal
nodes_areas_f <- drugs_N_filtered %>%
  distinct(therapeutic_area) %>%
  rename(name = therapeutic_area) %>%
  mutate(type = "Terápiás terület", atc_subgroup = NA_character_)

nodes_substances_f <- drugs_N_filtered %>%
  distinct(active_substance) %>%
  rename(name = active_substance) %>%
  mutate(type = "Hatóanyag") %>%
  left_join(atc_subgroups, by = c("name" = "active_substance"))

nodes_f <- bind_rows(nodes_areas_f, nodes_substances_f)

# Hálózat létrehozása súlyozott élekkel
network_f <- tbl_graph(nodes = nodes_f, edges = edges_weighted, directed = FALSE)

# Fokszám és klaszter hozzáadása
network_f <- network_f %>%
  activate(nodes) %>%
  mutate(
    degree = centrality_degree(),
    cluster = as.factor(group_louvain()))

Hálózat vizualizáció

# Hálózat vizualizációja - Interaktív ábra (jobb átláthatóság)
layout_coords <- create_layout(network_f, layout = "fr")
layout_coords$x <- layout_coords$x * 2
layout_coords$y <- layout_coords$y * 2

p <- ggraph(layout_coords) +
  geom_edge_link(aes(width = weight), alpha = 0.3, color = "gray50") +
  geom_point_interactive(aes(
    x = x, y = y,
    size = degree,
    color = ifelse(type == "Terápiás terület", "Terápiás terület", atc_subgroup),
    alpha = degree,
    tooltip = paste0(name, "\nKapcsolatok: ", degree),
    data_id = name
  )) +
  scale_size_continuous(range = c(2, 12)) +
  scale_alpha_continuous(range = c(0.5, 1)) +
  scale_edge_width_continuous(range = c(0.3, 3)) +
  scale_color_manual(
    values = c(
      "Terápiás terület" = "steelblue",
      "N01" = "#E41A1C",
      "N02" = "#FF7F00",
      "N03" = "#4DAF4A",
      "N04" = "#984EA3",
      "N05" = "#F781BF",
      "N06" = "#A65628",
      "N07" = "#999999"
    ),
    na.value = "coral"
  ) +
  labs(
    title = "Terápiás területek és hatóanyagok hálózata (Interaktív)",
    subtitle = "N (központi idegrendszer) ATC csoport - Kérem, vigye az egeret a pontokra a részletekért",
    color = "Típus / ATC alcsoport",
    size = "Kapcsolatok száma",
    edge_width = "Közös gyógyszerek"
  ) +
  theme_graph() +
  theme(legend.position = "right") +
  guides(alpha = "none")

girafe(
  ggobj = p,
  width_svg = 14,
  height_svg = 10,
  options = list(
    opts_zoom(max = 5),
    opts_hover(css = "fill:red;stroke:black;stroke-width:2;"),
    opts_tooltip(css = "background-color:white;padding:8px;border-radius:5px;border:1px solid gray;font-size:12px;")))
# Alapvető hálózati statisztikák
network_stats <- tibble(
  Mutató = c(
    "Csúcsok száma (nodes)",
    "Élek száma (edges)", 
    "Sűrűség (density)",
    "Átlagos fokszám (mean degree)",
    "Átmérő (diameter)",
    "Átlagos úthossz (mean distance)",
    "Tranzitivitás (clustering coefficient)"
  ),
  Érték = c(
    vcount(network_f),
    ecount(network_f),
    round(edge_density(network_f), 3),
    round(mean(degree(network_f)), 2),
    diameter(network_f),
    round(mean_distance(network_f), 2),
    round(transitivity(network_f), 3)))
network_stats
## # A tibble: 7 × 2
##   Mutató                                   Érték
##   <chr>                                    <dbl>
## 1 Csúcsok száma (nodes)                  110    
## 2 Élek száma (edges)                     108    
## 3 Sűrűség (density)                        0.018
## 4 Átlagos fokszám (mean degree)            1.96 
## 5 Átmérő (diameter)                       32    
## 6 Átlagos úthossz (mean distance)          6.61 
## 7 Tranzitivitás (clustering coefficient)   0
# Centralitás mutatók hozzáadása a hálózathoz
network_f <- network_f %>%
  activate(nodes) %>%
  mutate(
    degree = centrality_degree(),
    betweenness = round(centrality_betweenness(), 2),
    closeness = round(centrality_closeness(), 4),
    eigenvector = round(centrality_eigen(), 3))

# Top 20 csúcs különböző centralitás mutatók alapján
network_f %>%
  activate(nodes) %>%
  as_tibble() %>%
  arrange(desc(degree)) %>%
  select(name, type, degree, betweenness, closeness, eigenvector) %>%
  head(20)
## # A tibble: 20 × 6
##    name                           type  degree betweenness closeness eigenvector
##    <chr>                          <chr>  <dbl>       <dbl>     <dbl>       <dbl>
##  1 Parkinson Disease              Terá…     14       167.     0.0345       1    
##  2 amifampridine phosphate        Ható…     13       150      0.0312       0    
##  3 Schizophrenia                  Terá…     11        69.5    0.0625       0    
##  4 Epilepsy                       Terá…     10       176.     0.02         0    
##  5 Migraine Disorders             Terá…      6        15      0.167        0    
##  6 Alzheimer Disease              Terá…      5        54.3    0.0204       0.234
##  7 Opioid-Related Disorders       Terá…      5        10      0.2          0    
##  8 Pain                           Terá…      5        51.5    0.0172       0    
##  9 duloxetine                     Ható…      5        93      0.0145       0    
## 10 Bipolar Disorder               Terá…      4        13.5    0.0333       0    
## 11 Epilepsies, Partial            Terá…      4        41.5    0.0125       0    
## 12 sodium oxybate                 Ható…      4        12      0.125        0    
## 13 Sleep Initiation and Maintena… Terá…      3         5      0.2          0    
## 14 Depressive Disorder, Major     Terá…      3        41      0.0116       0    
## 15 Cancer                         Terá…      3        78.5    0.0263       0    
## 16 Neuralgia                      Terá…      3        63      0.0167       0    
## 17 Narcolepsy                     Terá…      3         9      0.111        0    
## 18 Pain, Postoperative            Terá…      3         3      0.333        0    
## 19 rivastigmine                   Ható…      3        37      0.0263       0.365
## 20 pregabalin                     Ható…      3       118      0.0189       0

Hálózatemezés mutatók

A kétoldalú (bipartite) hálózat 110 csúcspontot (terápiás területek és hatóanyagok) és 108 élt (kapcsolatot) tartalmaz. A hálózat ritkának tekinthető (sűrűség = 0.018), ami azt jelenti, hogy az összes lehetséges kapcsolat mindössze 1.8%-a realizálódik. Az átlagos fokszám 1.96, ami arra utal, hogy egy csúcspont átlagosan közel két másik csúcsponttal áll kapcsolatban.

A hálózat átmérője 32, az átlagos úthossz pedig 6.61, ami viszonylag szétszórt, kevéssé integrált struktúrára utal. A tranzitivitás (klaszterezési együttható) értéke 0, ami a bipartite hálózatok sajátossága, ahol azonos típusú csúcspontok nem kapcsolódnak közvetlenül egymáshoz.

Centralitás elemzés

A fokszám-centralitás alapján a legkiemelkedőbb terápiás területek a Parkinson-kór (fokszám = 14), a Szkizofrénia (fokszám = 11) és az Epilepszia (fokszám = 10). Ezek a betegségterületek számos különböző hatóanyaggal állnak kapcsolatban, ami a terápiás lehetőségek sokféleségét tükrözi.

A közöttiség-centralitás (betweenness) tekintetében az Epilepszia (176.0) és a Parkinson-kór (167.0) kiemelkedik, ami arra utal, hogy ezek a terápiás területek híd szerepet töltenek be a hálózatban, összekapcsolva más, egyébként elszigetelt klasztereket.

A hatóanyagok közül az amifampridine-foszfát rendelkezik a legmagasabb fokszámmal (13), amely a Lambert-Eaton myastheniás szindróma és más neuromuszkuláris betegségek kezelésében használatos. A duloxetin (fokszám = 5) és a pregabalin (fokszám = 3) szintén kiemelkedő pozíciót foglal el, ami széles terápiás alkalmazhatóságukat tükrözi (fájdalom, depresszió, szorongás, neuropátia).

Regressziós modellek

# Regressziós modell előkészítése
model_df <- top4_atc %>%
  transmute(
    revision_number,
    days_on_market,
    atc_group,
    generic,
    orphan_medicine,
    accelerated_assessment
  ) %>%
  filter(
    !is.na(revision_number),
    !is.na(days_on_market),
    !is.na(atc_group)
  ) %>%
  mutate(
    atc_group = factor(atc_group),
    generic = factor(generic),
    orphan_medicine = factor(orphan_medicine),
    accelerated_assessment = factor(accelerated_assessment),
    days_on_market_z = as.numeric(scale(days_on_market)))

skimr::skim(model_df)
Data summary
Name model_df
Number of rows 1026
Number of columns 7
_______________________
Column type frequency:
factor 4
numeric 3
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
atc_group 0 1 FALSE 4 L: 418, J: 245, A: 190, N: 173
generic 0 1 FALSE 2 FAL: 837, TRU: 189
orphan_medicine 0 1 FALSE 2 FAL: 931, TRU: 95
accelerated_assessment 0 1 FALSE 2 FAL: 996, TRU: 30

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
revision_number 0 1 15.62 12.45 0.00 6.00 13.00 22.00 89.00 ▇▃▁▁▁
days_on_market 0 1 3886.63 2542.55 124.00 1816.50 3462.50 5488.25 9968.00 ▇▇▆▂▃
days_on_market_z 0 1 0.00 1.00 -1.48 -0.81 -0.17 0.63 2.39 ▇▇▆▂▃

A revíziók számát befolyásoló tényezők azonosítására többváltozós lineáris regressziós modelleket építettem. A függő változó a revíziók száma volt, a prediktorok pedig az ATC csoport, a piacon töltött idő (z-standardizált), a generikus státusz, az orphan gyógyszer státusz és a gyorsított értékelés.

Az elemzés 1026 megfigyelésen alapult a hiányzó értékek kizárása után. A piacon töltött idő változót z-standardizáltam az együtthatók értelmezhetősége érdekében.

Két modellt hasonlítottam össze: hagyományos legkisebb négyzetek (OLS) regressziót és robusztus regressziót (MM-becslő), amely kevésbé érzékeny a kiugró értékekre.

Regressziós modellek létrehozása

# Regressziós modellek létrehozása - két modell: lm vs. lmrob (robusztus) 
mod_lm <- lm(
  revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment,
  data = model_df)

mod_rob <- lmrob(
  revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment,
  data = model_df)

# Regressziós modellek eredményei
summary(mod_lm)
## 
## Call:
## lm(formula = revision_number ~ atc_group + days_on_market_z + 
##     generic + orphan_medicine + accelerated_assessment, data = model_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.208  -4.405  -0.668   4.441  62.413 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 13.7840     0.7468  18.458  < 2e-16 ***
## atc_groupJ                   3.1040     0.9543   3.253  0.00118 ** 
## atc_groupL                   3.6263     0.8707   4.165 3.38e-05 ***
## atc_groupN                   2.2845     1.0407   2.195  0.02838 *  
## days_on_market_z             7.2030     0.3300  21.825  < 2e-16 ***
## genericTRUE                 -4.1782     0.8287  -5.042 5.46e-07 ***
## orphan_medicineTRUE         -1.0943     1.1333  -0.966  0.33449    
## accelerated_assessmentTRUE   3.6924     1.8587   1.986  0.04725 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.828 on 1018 degrees of freedom
## Multiple R-squared:  0.3808, Adjusted R-squared:  0.3765 
## F-statistic: 89.42 on 7 and 1018 DF,  p-value: < 2.2e-16
summary(mod_rob)
## 
## Call:
## lmrob(formula = revision_number ~ atc_group + days_on_market_z + generic + 
##     orphan_medicine + accelerated_assessment, data = model_df)
##  \--> method = "MM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.543  -3.935  -0.299   4.045  60.442 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 14.2598     0.4997  28.536  < 2e-16 ***
## atc_groupJ                   3.6720     0.8605   4.267 2.16e-05 ***
## atc_groupL                   2.9631     0.6137   4.828 1.59e-06 ***
## atc_groupN                   1.5387     0.6807   2.261  0.02399 *  
## days_on_market_z             8.8971     0.4217  21.099  < 2e-16 ***
## genericTRUE                 -3.5214     0.4465  -7.886 7.98e-15 ***
## orphan_medicineTRUE         -0.1621     0.5224  -0.310  0.75643    
## accelerated_assessmentTRUE   3.6881     1.1926   3.092  0.00204 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 6.065 
## Multiple R-squared:  0.6258, Adjusted R-squared:  0.6232 
## Convergence in 20 IRWLS iterations
## 
## Robustness weights: 
##  28 observations c(33,45,238,257,293,302,303,305,411,501,976,990,991,992,993,994,995,1014,1015,1018,1019,1020,1021,1022,1023,1024,1025,1026)
##   are outliers with |weight| = 0 ( < 9.7e-05); 
##  103 weights are ~= 1. The remaining 895 ones are summarized as
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.001371 0.828100 0.953900 0.850500 0.986900 0.999000 
## Algorithmic parameters: 
##        tuning.chi                bb        tuning.psi        refine.tol 
##         1.548e+00         5.000e-01         4.685e+00         1.000e-07 
##           rel.tol         scale.tol         solve.tol          zero.tol 
##         1.000e-07         1.000e-10         1.000e-07         1.000e-10 
##       eps.outlier             eps.x warn.limit.reject warn.limit.meanrw 
##         9.747e-05         4.351e-12         5.000e-01         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##            500             50              2              1            200 
##    maxit.scale      trace.lev            mts     compute.rd fast.s.large.n 
##            200              0           1000              0           2000 
##                   psi           subsampling                   cov 
##            "bisquare"         "nonsingular"         ".vcov.avar1" 
## compute.outlier.stats 
##                  "SM" 
## seed : int(0)

Regressziós modellek összehasolítása

# Együtthatók összehasonlítása
coef_comparison <- bind_rows(
  tidy(mod_lm, conf.int = TRUE) %>% mutate(model = "OLS"),
  tidy(mod_rob, conf.int = TRUE) %>% mutate(model = "Robust"))

coef_comparison %>%
  select(model, term, estimate, std.error, conf.low, conf.high, p.value) %>%
  mutate(
    estimate = round(estimate, 3),
    std.error = round(std.error, 3),
    conf.low = round(conf.low, 3),
    conf.high = round(conf.high, 3),
    p.value = format_p(p.value))
## # A tibble: 16 × 7
##    model  term                     estimate std.error conf.low conf.high p.value
##    <chr>  <chr>                       <dbl>     <dbl>    <dbl>     <dbl> <chr>  
##  1 OLS    (Intercept)                13.8       0.747   12.3      15.2   < 0.001
##  2 OLS    atc_groupJ                  3.10      0.954    1.23      4.98  0.001  
##  3 OLS    atc_groupL                  3.63      0.871    1.92      5.34  < 0.001
##  4 OLS    atc_groupN                  2.28      1.04     0.242     4.33  0.028  
##  5 OLS    days_on_market_z            7.20      0.33     6.56      7.85  < 0.001
##  6 OLS    genericTRUE                -4.18      0.829   -5.80     -2.55  < 0.001
##  7 OLS    orphan_medicineTRUE        -1.09      1.13    -3.32      1.13  0.334  
##  8 OLS    accelerated_assessmentT…    3.69      1.86     0.045     7.34  0.047  
##  9 Robust (Intercept)                14.3       0.5     13.3      15.2   < 0.001
## 10 Robust atc_groupJ                  3.67      0.86     1.98      5.36  < 0.001
## 11 Robust atc_groupL                  2.96      0.614    1.76      4.17  < 0.001
## 12 Robust atc_groupN                  1.54      0.681    0.205     2.87  0.024  
## 13 Robust days_on_market_z            8.90      0.422    8.07      9.72  < 0.001
## 14 Robust genericTRUE                -3.52      0.447   -4.40     -2.65  < 0.001
## 15 Robust orphan_medicineTRUE        -0.162     0.522   -1.19      0.862 0.756  
## 16 Robust accelerated_assessmentT…    3.69      1.19     1.35      6.03  0.002
# Együtthatók vizualizációja
ggplot(coef_comparison %>% filter(term != "(Intercept)"), 
       aes(x = estimate, y = term, color = model)) +
  geom_point(position = position_dodge(width = 0.5), size = 3) +
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), 
                 position = position_dodge(width = 0.5), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  labs(
    title = "Regressziós együtthatók összehasonlítása",
    subtitle = "OLS vs. Robusztus regresszió",
    x = "Becsült együttható (95% CI)",
    y = "Prediktor",
    color = "Modell"
  ) +
  theme_minimal()

Regressziós eredmények

OLS modell

Az OLS modell szignifikánsan illeszkedett az adatokra, F(7, 1018) = 89.42, p < .001, és a revíziószám varianciájának 38.1%-át magyarázta (R² = .381, R²adj = .377).

A piacon töltött idő volt a legerősebb prediktor: minden egységnyi növekedés a z-standardizált piacon töltött időben 7.20 egységgel növelte a revíziók számát, B = 7.20, 95% CI [6.56, 7.85], p < .001. Ez azt jelenti, hogy a régebben engedélyezett gyógyszerek több revízión estek át.

Az ATC csoport szignifikáns prediktornak bizonyult. A referenciakategóriához (A csoport) képest a J csoport (B = 3.10, p = .001), az L csoport (B = 3.63, p < .001) és az N csoport (B = 2.28, p = .028) egyaránt szignifikánsan magasabb revíziószámmal rendelkezett.

A generikus státusz szignifikáns negatív összefüggést mutatott a revíziók számával (B = -4.18, 95% CI [-5.80, -2.55], p < .001), ami arra utal, hogy a generikus gyógyszerek átlagosan 4.18-cal kevesebb revízión esnek át, mint az originális készítmények.

A gyorsított értékelés szignifikáns pozitív prediktora volt a revíziószámnak (B = 3.69, 95% CI [0.05, 7.34], p = .047). Az orphan gyógyszer státusz nem mutatott szignifikáns összefüggést (B = -1.09, p = .334).

Robusztus modell

A robusztus regresszió (MM-becslő) 28 megfigyelést azonosított kiugró értékként. A robusztus modell lényegesen magasabb magyarázóerővel rendelkezett (R² = .626, R²adj = .623), ami arra utal, hogy a kiugró értékek jelentősen torzították az OLS becsléseket.

A prediktorok mintázata hasonló maradt, de néhány fontos eltérés mutatkozott. A piacon töltött idő hatása erősebb volt a robusztus modellben (B = 8.90, 95% CI [8.07, 9.72], p < .001). Illetve a gyorsított értékelés hatása is robusztusabbá vált (B = 3.69, p = .002).

Modell összehasonlítás

A keresztvalidáció (80-20% train-test split) alapján a robusztus modell kissé jobb prediktív teljesítményt mutatott: RMSE = 11.5 vs. 11.7, MAE = 7.32 vs. 7.68, R²pred = .376 vs. .356.

A Cook-féle távolság alapján 51 megfigyelés (5.0%) minősült potenciális kiugró értéknek (Cook’s d > 4/n). A residuumok vizsgálata heteroszkedaszticitásra utalt, ami a robusztus módszerek alkalmazását indokolttá tette.

Regressziós modellek illeszkedése

# Modell illeszkedési mutatók
compare_performance(mod_lm, mod_rob, metrics = c("R2", "R2_adjusted", "RMSE", "AIC", "BIC"))
## # Comparison of Model Performance Indices
## 
## Name    | Model |    R2 |  RMSE |  AIC (weights) |  BIC (weights)
## -----------------------------------------------------------------
## mod_lm  |    lm | 0.381 | 9.789 | 7610.9 (>.999) | 7655.3 (>.999)
## mod_rob | lmrob | 0.626 | 9.956 |        (>.999) |        (>.999)

Regressziós modellek - Keresztvalidáció

# Regressziós modellek keresztvalidációja
set.seed(123)
idx <- sample(seq_len(nrow(model_df)), size = floor(0.8 * nrow(model_df)))
train <- model_df[idx, ]
test  <- model_df[-idx, ]

m1 <- lm(
  revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment, 
  data = train)
m2 <- lmrob(
  revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment, 
  data = train)

pred1 <- predict(m1, newdata = test)
pred2 <- predict(m2, newdata = test)

# Predikciós hibák
rmse <- function(y, yhat) sqrt(mean((y - yhat)^2))
mae  <- function(y, yhat) mean(abs(y - yhat))
r2_pred <- function(y, yhat) 1 - sum((y - yhat)^2) / sum((y - mean(y))^2)

prediction_comparison <- tibble(
  Modell = c("OLS", "Robust (lmrob)"),
  RMSE = round(c(rmse(test$revision_number, pred1), rmse(test$revision_number, pred2)), 3),
  MAE = round(c(mae(test$revision_number, pred1), mae(test$revision_number, pred2)), 3),
  R2_pred = round(c(r2_pred(test$revision_number, pred1), r2_pred(test$revision_number, pred2)), 3))
prediction_comparison
## # A tibble: 2 × 4
##   Modell          RMSE   MAE R2_pred
##   <chr>          <dbl> <dbl>   <dbl>
## 1 OLS             11.7  7.68   0.356
## 2 Robust (lmrob)  11.5  7.32   0.376

Regressziós modellek - Diagnosztika

# Regressziós modellek diagnosztikája (ábra)
check_model(mod_lm)

check_model(mod_rob)

Regressziós modellek - Kiugró értékek

# Cook's distance az OLS modellhez
model_df$cooks_d <- cooks.distance(mod_lm)
model_df$fitted <- fitted(mod_lm)
model_df$residuals <- residuals(mod_lm)

# Kiugró értékek azonosítása (Cook's d > 4/n)
n <- nrow(model_df)
outliers <- model_df %>%
  filter(cooks_d > 4/n)

cat("Kiugró értékek száma:", nrow(outliers), "\n")
## Kiugró értékek száma: 51
cat("Kiugró értékek aránya:", round(nrow(outliers)/n * 100, 1), "%\n")
## Kiugró értékek aránya: 5 %
# Vizualizáció
ggplot(model_df, aes(x = fitted, y = residuals)) +
  geom_point(aes(size = cooks_d, color = cooks_d > 4/n), alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_color_manual(values = c("FALSE" = "gray50", "TRUE" = "red")) +
  labs(
    title = "Residuumok vs. Illesztett értékek",
    subtitle = "Piros pontok: kiugró értékek (Cook's d > 4/n)",
    x = "Illesztett értékek",
    y = "Residuumok",
    size = "Cook's d",
    color = "Kiugró"
  ) +
  theme_minimal()

Regressziós modellek - Összefoglaló táblázat

# Összefoglaló táblázat (modelsummary)

modelsummary(
  list("OLS" = mod_lm, "Robust" = mod_rob),
  estimate = "{estimate} ({std.error})",
  statistic = "p.value",
  stars = TRUE,
  title = "Regressziós modellek összehasonlítása",
  coef_rename = c(
    "atc_groupJ" = "ATC: J",
    "atc_groupL" = "ATC: L",
    "atc_groupN" = "ATC: N",
    "days_on_market_z" = "Piacon töltött idő (z-score)",
    "genericTRUE" = "Generikus",
    "orphan_medicineTRUE" = "Orphan gyógyszer",
    "accelerated_assessmentTRUE" = "Gyorsított eljárás"
  )
)
Regressziós modellek összehasonlítása
OLS Robust
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
(Intercept) 13.784 (0.747) 14.260 (0.500)
(<0.001) (<0.001)
ATC: J 3.104 (0.954) 3.672 (0.860)
(0.001) (<0.001)
ATC: L 3.626 (0.871) 2.963 (0.614)
(<0.001) (<0.001)
ATC: N 2.285 (1.041) 1.539 (0.681)
(0.028) (0.024)
Piacon töltött idő (z-score) 7.203 (0.330) 8.897 (0.422)
(<0.001) (<0.001)
Generikus -4.178 (0.829) -3.521 (0.447)
(<0.001) (<0.001)
Orphan gyógyszer -1.094 (1.133) -0.162 (0.522)
(0.334) (0.756)
Gyorsított eljárás 3.692 (1.859) 3.688 (1.193)
(0.047) (0.002)
Num.Obs. 1026 1026
R2 0.381 0.626
R2 Adj. 0.377 0.623
AIC 7610.9
BIC 7655.3
Log.Lik. -3796.426
F 89.422
RMSE 9.79 9.96

Összegzés és konklúzió

Jelen elemzés az Európai Gyógyszerügynökség gyógyszer-engedélyezési adatbázisát vizsgálta négy fő kutatási kérdés mentén.

Engedélyezési trendek. A gyógyszer-engedélyezések száma jelentős fluktuációt mutatott 1995 és 2023 között, kiugró növekedéssel a COVID-19 pandémia időszakában. Az L csoport (daganatellenes szerek) dominanciája különösen az elmúlt évtizedben vált hangsúlyossá.

ATC csoportok közötti különbségek. Szignifikáns különbségek mutatkoztak a revíziók számában (H(3) = 9.39, p = .025) és a piacon töltött időben (H(3) = 48.2, p < .001) a négy leggyakoribb ATC csoport között. Az L csoport gyógyszerei szignifikánsan rövidebb ideje vannak a piacon, ami feltételezhetően az onkológiai területen tapasztalható gyorsabb innovációval magyarázható.

Revíziószámot befolyásoló tényezők. A regressziós elemzés alapján a piacon töltött idő, az ATC csoport, a generikus státusz és a gyorsított értékelés szignifikáns prediktorai a revíziók számának. A robusztus regresszió magasabb magyarázóerőt mutatott (R² = .63 vs. .38), ami a kiugró értékek jelentős torzító hatására utal.

Hálózati kapcsolatok. A központi idegrendszeri gyógyszerek hálózatelemzése feltárta a terápiás területek és hatóanyagok közötti komplex kapcsolatrendszert. A Parkinson-kór, a Szkizofrénia és az Epilepszia központi pozíciót foglal el a hálózatban, ami ezeknek a betegségterületeknek a terápiás diverzitását tükrözi.

Limitációk

Az elemzés korlátai közé tartozik, hogy az adatbázis csak az EMA által centralizált eljárásban engedélyezett gyógyszereket tartalmazza, így az eredmények nem általánosíthatók a nemzeti szinten engedélyezett készítményekre. A revíziók száma mint kimeneti változó korlátozottan tükrözi a gyógyszerek tényleges biztonságossági és hatékonysági profilját. A hálózatelemzés keresztmetszeti jellegű, és nem tárja fel az időbeli változásokat a kapcsolatrendszerben.

Jövőbeli irányok

További kutatások vizsgálhatnák a revíziók tartalmát és típusát, a forgalomból történő kivonások prediktorait, valamint az engedélyezési folyamat időtartamának meghatározó tényezőit.